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The spectrum of ideal magnetohydrodynamic (MHD) 
pressure-driven (ballooning) modes in strongly nonaxisym- 
metric toroidal systems is difficult to analyze numerically ow- 
ing to the singular nature of ideal MHD caused by lack of 
an inherent scale length. In this paper, ideal MHD is regu- 
larized by using a fc-space cutoff, making the ray tracing for 
the WKB ballooning formalism a chaotic Hamiltonian bil- 
liard problem. The minimum width of the toroidal Fourier 
spectrum needed for resolving toroidally localized ballooning 
modes with a global eigenvalue code is estimated from the 
Weyl formula. This phase-space-volume estimation method 
is applied to two stellarator cases. 

PACS numbers: 52.35. Py, 52.55.Hc, 05.45. Mt 

In design studies for new magnetic confinement devices 
for fusion plasma experiments (e.g. investigations [QJ|] 
leading to the proposed National Compact Stellarator 
Experiment, NCSX J3|), the maximum pressure that can 
stably be confined in any proposed magnetic field con- 
figuration is routinely estimated by treating the plasma 
as an ideal magnetohydrodynamic (MHD) fluid. One 
linearizes about a sequence of equilibrium states with in- 
creasing pressure, and studies the spectrum of normal 
modes (frequency ui) to determine when there is a com- 
ponent with Imw > 0, signifying instability. 

Even with the simplification obtained by using the 
ideal MHD model, the computational task of determining 
the theoretical stability of a three-dimensional (i.e. non- 
axisymmetric) device, such as NCSX or the four currently 
operating helical axis stellators Q , remains a challenging 
one. 

The problem can be posed as a Lagrangian field the- 
ory, with the potential term being the energy functional 
SW ||. For a static equilibrium, the kinetic energy is 
quadratic in w, so that uj 2 is real. Thus instability oc- 
curs when ui 2 < 0. There are two main approaches to 
analyzing the spectrum — local and global. 
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In the local approach, which is used for analytical sim- 
plification, one orders the scale length of variation of the 
eigenfunction across the magnetic field lines to be short 
compared with equilibrium scale lengths Q . Both inter- 
change and ballooning stability can be treated by solving 
the general ballooning equations , a system of ordinary 
differential equations defined on a given magnetic field 
line. 

The global (Galerkin) approach is to expand the 
plasma displacement field in a finite basis set, inserting 
this ansatz in the Lagrangian to find a matrix eigen- 
value representation of the spectral problem. This ap- 
proach has been implemented for ideal MHD in three- 
dimensional plasmas in two codes, TERPSICHORE ||] 
and CAS3D §. 

Although the Galerkin approach is potentially exact, if 
one could use a complete, infinite basis set, it is in prac- 
tice computationally challenging due to the large number 
of basis functions required to resolve localized instabili- 
ties. This leads to very large matrices which must be 
diagonalized by iterative methods. There is a need for 
analytical insight to determine a suitable truncated ba- 
sis set and to predict the nature of the spectrum, e.g. 
whether it is continuous or discrete. 

Such insight may be obtained by a hybrid local- 
global approach, in which one uses a Wentzel-Kramers- 
Brillouin (WKB) representation of the eigenfunction. In 
the short-wavelength limit, the same analytical simplifi- 
cations as are obtained in the local approach are found to 
give a local dispersion relation that can be used to give 
information on the global spectrum by using ray tracing 
and semiclassical quantization. 

In axisymmetric systems [fiof or in cases where helical 
ripple can be averaged out, giving an adiabatic invariant, 
pl|,^2[ , the ray equations are integrable and hence the 
spectrum is characterized by "good quantum numbers" . 

However, it has been known for many years |?J that the 
ray-tracing problem in strongly three-dimensional sys- 
tems is singular because, in the absence of an adiabatic 
invariant, the phase-space motion is not bounded — the 
rays escape to infinity in the wavevector sector. Dewar 
and Glasser Q argued that this gives rise to a contin- 
uous unstable spectrum, with correspondingly singular 
generalized eigenfunctions. (A more rigorous treatment 
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involves the concept of the essential spectrum and Weyl 
sequences 

Our proposed regularization of this singularity can be 
understood using a simple quantum analogy. Consider 
the one-dimensional time-independent Schrodinger equa- 
tion Hip = Eip in the limit as the mass of the par- 
ticle goes to infinity. Then the kinetic energy disap- 
pears and the Hamiltonian becomes H = V(x), where 
V is the potential energy, assumed here to be the har- 
monic oscillator potential, \x 2 in suitable units. In the 
usual Hilbert space the energy spectrum is continuous: 
E > and the (generalized) eigenfunctions singular: 
ip(x) = S(x — xe) ± 8(x + xe), where V(xe) = E. 

We now seek a regularization of this problem by re- 
stricting ip to the space of functions with a finite band- 
width in wavenumber k: 
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This truncated Fourier-integral representation models 
what occurs when one seeks to find the spectrum numer- 
ically using a truncated Fourier-series representation. 

We take as starting point a Lagrangian for the wave- 
function, 



L 



ip*[E -V(x)]4>dx 
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Inserting Eq. (|l|) in Eq. 
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This is infinite unless we require the coefficients of the 
(5-functions to vanish. That is, ip k = at k = ±fc max . 
The Euler-Lagrange equation is (d 2 /2dk 2 + E)ip k = 0, 
which has the solutions exp ±i(2E) 1 ' 2 k, These waves 
would propagate to infinity if it were not for the reflecting 
boundary conditions at ±fc max we have just derived. 

That is, we have removed the continuum by box quan- 
tization in fc-space. In the following we shall do the same 
for the ballooning mode problem. 

As in we write the magnetic field of an arbi- 
trary three-dimensional toroidal equilibrium plasma with 
nested magnetic flux surfaces labeled by an arbitrary pa- 
rameter s as B = VCx Vip - qVOxVip = VaxVV>, 
where a = £ — qd. Here, 9 and C are the poloidal and 
toroidal angles, respectively, ip(s) is the poloidal flux 
function, and q(s) is the inverse of the rotational trans- 
form. Since B-Vs = B-Va = 0, s and a serve to label 
an individual field line. 

We take the stream function H to be given by ip = 
<pexp(iS — iojt) , where <p(0\s, a) is assumed to vary on 
the equilibrium scale. The phase variation is taken to be 



rapid, so k = V5 is ordered to be large. The frequency 
ui is ordered 0(1), which requires that the wave vector be 
perpendicular to B: k-B = 0. (In this study we consider 
unstable ideal MHD modes, ui 2 < 0.) 

It immediately follows that the eikonal is constant on 
each field line: S = S(a,s). From the definition of the 
wave vector, k = k a Va + fc s Vs = k a \Va + 6'fe(7 / (s)Vs] 
where k a = dS/da and k s = dS/ds. Here the anglelike 
ballooning parameter 9 k appears naturally as the ratio 
k s /q'(s)k a @. 

The ballooning equation emerges in the large |k| ex- 
pansion |7],|6| as an ordinary differential equation to be 
solved on each field line {a, s) with given (k a ,k s ) under 
the boundary condition @(0) — ► at infinity to give the 
eigenvalue A(a, s, k a , k s ). This constitutes a local disper- 
sion relation A = puj 2 (the mass density p being assumed 
constant everywhere). 

The ray equations are the characteristics of the eikonal 
equation X(a, s, d a S, d s S) = pui 2 . These are Hamiltonian 
equations of motion with a, s the generalized coordinates, 
k a , k s the canonically conjugate momenta, and A as the 
Hamiltonian. 

In axi- or helically symmetric systems all field lines on 
a given magnetic surface are equivalent — a is ignorable 
and k a is a constant of the motion. In this case the equa- 
tions are integrable and semiclassical quantization can be 
used to predict the approximate spectrum of global bal- 
looning instabilities [ jlOj ] . This technique can sometimes 
be applied successfully, even in nonsymmetric systems, if 
there are regions of phase space with a large measure of 
invariant tori JT^,[TT|] . In [jllj this was verified using the 
global eigenvalue code TERPSICHORE @. 

At the other extreme, if the ray orbits are chaotic (but 
still bounded) then the global spectrum is not regularly 
structured, but must rather be described statistically by 
the density of states and the probability distribution of 
level spacings using the techniques of quantum chaos the- 
ory (see e.g. |l6|Jl7|1 ). 

However, because of the scale invariance of the ideal 
MHD equations, A depends only on the direction of k, 
not on its magnitude: A = A(a, s, 6 k )- This has the 
consequence that the ray orbits are unbounded in phase 
space, so, strictly speaking, ideal MHD gives rise to a 
quantum chaotic scattering [ ^6|Jr^| problem rather than 
a straight quantum chaos problem. This leads to the con- 
tinuous spectrum |7j] with singular generalized eigenfunc- 
tions that cannot really be represented using the simple 
eikonal ansatz. 

On the other hand, the absence of a natural length 
scale in ideal MHD is a mathematical artifact. Phys- 
ically, the ion Larmor radius provides a lower cutoff in 
space, or an upper cutoff in |k|, beyond which ideal MHD 
ceases to apply. The ballooning equation is also physi- 
cally regularized by inclusion of diamagnetic drift |l8|,|l5| . 

However, since in general it leads to a complex ray 
tracing problem Q , we shall not attempt to model dia- 
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magnetic drift stabilization in this paper. Rather, we 
regularize the ray equations simply by adding a barrier 
term to the effective ray "Hamiltonian" H(a, s,k a ,k s ), 



H = A(a, s, k a , k s ) + U{k a ) 



(4) 



where the barrier potential we use is U(k a ) = K(\k a \ — 
W) 2 for | k a | > /c max and for \k a \ < fc max . In the limit 
of the constant K — > oo, this infinite box potential gives 
the ideal MHD ray equations for I k a \ < fc max and reflect- 



ing boundary conditions at \k D 



Thus we have a 



two-degree of freedom Hamiltonian billiard problem. 

Although overly crude for modeling FLR stabilization, 
the cutoff at \k a \ = fc max provides a reasonable model for 
representing the finite spectral bandwidth in the toroidal 
Fourier mode number (n) representation used in the 
global eigenvalue codes TERPSICHORE §] and CAS3D 
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FIG. 1. The sections 6 k = and q = 0.893 of the topo- 
logically spherical isosurfaces of the central, (0,0), ballooning 
mode branch, bounded by the isosurface A = —6 (arbitrary 
units). The darker shades denote higher growth rates, the 
peak corresponding to A « —8. 

Using ballooning- unstable plasma equilibria calculated 
for the H-1NF heliac |(],f| using the VMEC code [fH, 
detailed parameter scans have been undertaken for two 
cases. The first case studied J|^] was obtained by increas- 
ing the pressure gradient of a marginally stable equilib- 
rium [p3| uniformly across the plasma and thus was bal- 
looning unstable at the edge of the plasma. The ray 
tracing problem for this case would involve consideration 
of the effect of the plasma boundary. 

Thus a second equilibrium, ballooning stable near the 
edge of the plasma, was calculated for the purposes of 
the present paper. This case has a more peaked pressure 
profile than the first, but both have average ~ 1%, 
where /3 is the ratio of plasma pressure to magnetic field 
pressure. 

The g-profiles are not monotonic — in the peaked pres- 
sure profile case studied in this paper, q was 0.8895 on 
the magnetic axis, rising to a maximum value of 0.8964 
quite close to the magnetic axis, then falling monoton- 
ically to 0.8675. Clearly the (global) magnetic shear is 
very weak. Despite this fact and the non-monotonicity, 
there is some formal simplification in choosing s = q, and 



we have taken s = q since the region of plasma studied 
is in a monotonic-decreasing part of the q-profile (the 
decreasing region outside the maximum-q surface). 

In these scans the most unstable ballooning eigen- 
value was tabulated on a three-dimensional grid in s, a, 9k 
space. The dependence on a was found to be rapid. The 
dependence on 9k was much slower, but the variation was 
sufficient that the higher-growth-rate isosurfaces formed 
a set of distinct, topologically spherical branches. It was 
argued in |2^] that this branch structure is produced by 
Anderson localization in bad curvature regions due to 
the strong breaking of both helical and axisymmetry in 
H-1NF. 

According to the perturbation expansion in q' de- 
scribed in j22|, a quadratic form in a, 9k should form a 
good approximation to A — X m i n (q) in the neighborhood 
of the central branch. Accordingly a least-squares fit on 
each surface was performed to provide a simple analytical 
description of the (0,0) branch. 

The radial dependence of the fitting coefficients was 
approximated by fitting to third-degree polynomials in 
q. Sections of the resulting approximation to the cen- 
tral branch are shown in Fig. |l[ The isosurface spans a 
substantial range of magnetic surfaces within the plasma 
- the narrow range of variation in q is due to the low 
magnetic shear in H-1NF. 

In order to establish the nature of the ray dynamics de- 
scribed by the regularized Hamiltonian, Eq. (^), a numer- 
ical integration with cutoff at fc max = 50 was performed 
with initial conditions q = q2, a — 0, and k a = 5, where 
[9i><72] = [0.8852,0.8951] is the g-range spanned by the 
A = —6 isosurface as seen in Fig. |l|. (A run with k a = 10 
was also performed, with similar results.) Choosing the 
value K = 1 gave a good compromise between the sharp 
boundary potential to be modeled, and the smooth po- 
tential required for the numerical integration. The orbit 
remained on the "energy shell" A = — 6 to within an ac- 
curacy of one part in 10 6 over the "time" interval of the 
integration, 7500. 
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FIG. 2. Two views of intersections with the Poincare sur- 
face of section a = 0. 

The two Poincare plots in Fig. show the orbit to be 
strongly chaotic, filling the "energy shell" ergodically, ex- 
cept that the regions k a > and k a < are dynamically 
disjoint. The solid curve shown surrounding the outer 
limits of the "energetically accessible" region is calcu- 
lated by solving A(0, q, A: 9 /A: max ) = —6. 
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According to the Weyl formula Q pp. 257-261], 
the number, iV(A max ), of global eigenmodes with eigen- 
values below the eigenvalue A max is given, asymp- 
totically in the limit N — > oo, as A^(A max ) ~ 
W4D(A ma x)/(27r) 2 . Here ^4D(A m ax) is the volume of the 
dynamically acessible 4-dimensional phase-space region 
\(a,q,k g /k a ) < A max , < k a < fc max . The k a integra- 
tion can be performed analytically, giving U4D(A ma x) = 
^max^DtAmax), where u 3 D(A max ) is the volume within 
the isosurface A(a, = A max . Thus 

N(X max ) ~ ^2^ ax t;3D(A max ) . (5) 

We can make a rather rough estimate of the minimum 
value of n max required for CAS3D or TERPSICHORE 
to find even one eigenvalue with A < A max by setting 
A(A max ) = 1 and calculating fc max w n max from Eq. (§). 
This gives n m ^(N = 1) ~ (87t 2 /w 3 d) 1/2 . 

The isosurface A = —6 studied above is about the 
largest of the disjoint topologically spherical isosurfaces 
corresponding to the highly toroidally localized strongly 
ballooning unstable regions of a, q, Ok space. (For A > —6 
the isosurfaces are no longer topologically spherical.) Us- 
ing the polynomial fits described above, we calculate 
w 3D (-6) = 0.02158. This gives n max (N = 1) « 60. As- 
suming that the dominant contributions to the MHD en- 
ergy 5W come from the rational surfaces intersecting the 
A = — 6 isosurface, we thus predict that it would be neces- 
sary to include, as a minimum set, basis functions corre- 
sponding to one of the two "mode families" II contained 
in the set (n,m) = (9,8), (18,16), (19,17), (27,24), 
(28,25), (35,31), (36,32), (37,33), (38,34), (44,39), 
(45,40), (46,41), (47,42), (53,47), (54,48), (55,49,), 
(56, 50), and (57, 51) to resolve a toroidally localized bal- 
looning mode. (Here n, m are the toroidal and poloidal 
Fourier mode numbers, respectively.) 

The large value of n max (iV = 1) required, and the un- 
usual spread in n required in the basis set, will make 
these modes difficult to resolve using global eigenvalue 
codes (e.g. the simplifying phase factor method some- 
times used in CAS3D studies M would not be appropri- 
ate). It is hoped that the Weyl formula estimate above 
will act as a guide in a future more extensive study us- 
ing such a code. Physically, the large value of n max sug- 
gests that toroidally localized ballooning modes in H-1NF 
should be subject to strong FLR stabilization. 

We can also apply the same approach to the toroidally 
localized ballooning branches found in the Large Heli- 
cal Device (LHD) study fil| . From the plots in |1^] we 
estimate u 3 d ~ 0.05, which gives n max (N = 1) w 40. 

The ballooning calculations were carried out on the 
Australian National University Supercomputer Facility's 
Fujitsu VPP300 vector processor. We thank Dr. H. J. 
Gardner for providing the H-l heliac VMEC input files 
and Dr. S. P. Hirshman for use of the VMEC equilibrium 
code. Some of this work was done while one of us (RLD) 



was a visiting scientist at Princeton University Plasma 
Physics Laboratory, supported under US DOE contract 
No. DE-AC02-76CH0-3703. Useful conversations with 
Drs. M. Redi and A.H. Boozer arc gratefully acknowl- 
edged. 



[1] A. H. Reiman et ai, Plasma Physics Reports 23, 472 
(1997). 

[2] A. Reiman et ai, Plasma Phys. Control. Fusion 41, B273 

(1999) . 

[3] G. H. Nielson et ai, Phys. Plasmas 7, 1911 (2000). 

[4] B. D. Blackwell, Bull. Am. Phys. Soc. 45, 289 (2000) 

[invited paper, to be published in Phys. Plasmas (2001)]. 
[5] I. B. Bernstein, E. A. Frieman, M. D. Kruskal, and R. M. 

Kulsrud, Proc. R. Soc. London Ser. A 244, 17 (1958). 
[6] R. L. Dewar, J. Plasma and Fusion Res. 73, 1123 (1997). 
[7] R. L. Dewar and A. H. Glasser, Phys. Fluids 26, 3038 

(1983). 

[8] D. V. Anderson et ai, Int. J. Supercomp. Appl. 4, 34 
(1990). 

[9] C. Schwab, Phys. Fluids B 5, 3195 (1993). 
[10] R. L. Dewar, J. Manickam, R. C. Grimm, and M. S. 

Chance, Nucl. Fusion 21, 493 (1981), corrigendum: Nucl. 

Fusion, 22 (1982) 307. 
[11] W. A. Cooper, D. B. Singleton, and R. L. Dewar, Phys. 

Plasmas 3, 275 (1996), erratum: Phys. Plasmas 3, 3520 

(1996). 

[12] P. Cuthbert et ai, Phys. Plasmas 5, 2921 (1998). 

[13] E. Hameiri, Commun. Pure Appl. Math. 38, 43 (1985). 

[14] A. E. Lifschitz, Magnetohydrodynamics and Spectral The- 
ory (Kluwer, Dordrecht, The Netherlands, 1989), pp. 
416-423. 

[15] W. M. Nevins and L. D. Pearlstein, Phys. Fluids 31, 1988 
(1988). 

[16] M. C. Gutzwiller, Chaos in Classical and Quantum Me- 
chanics, Interdisciplinary Applied Mathematics Series, 
Vol. 1 (Springer-Verlag, New York, 1990). 

[17] E. Ott, Chaos in Dynamical Systems (Cambridge Univ. 
Press, Cambridge, U.K., 1993). 

[18] W. M. Tang, R. L. Dewar, and J. Manickam, Nucl. Fusion 
22, 1079 (1982). 

[19] R. J. Hastie, P. J. Catto, and J. J. Ramos, Bull. Am. 
Phys. Soc. 45, 363 (2000). 

[20] S. M. Hamberger, B. D. Blackwell, L. E. Sharp, and D. B. 
Shenton, Fusion Technol. 17, 123 (1990). 

[21] S. P. Hirshman and O. Betancourt, J. Comput. Phys. 96, 
99 (1991). 

[22] P. Cuthbert and R. L. Dewar, Phys. Plasmas 7, 2302 

(2000) . 

[23] W. A. Cooper and H. J. Gardner, Nucl. Fusion 34, 729 
(1994). 



4 



